{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example creating the stiffness for a linear elasticity element using tensors or matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using JuAFEM\n",
    "using Tensors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stiffness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stiffness using normal matrices\n",
    "function ke_element_mat!{T, dim}(Ke, X::Vector{Vec{dim, T}}, fe_values::CellScalarValues{dim}, Ee, B, DB, BDB)\n",
    "    n_basefuncs = getnbasefunctions(fe_values)\n",
    "    @assert length(X) == n_basefuncs\n",
    "    \n",
    "    reinit!(fe_values, X)\n",
    "    for q_point in 1:getnquadpoints(fe_values)\n",
    "        for i in 1:n_basefuncs\n",
    "            dNdx = shape_gradient(fe_values, q_point, i)[1]\n",
    "            dNdy = shape_gradient(fe_values, q_point, i)[2]\n",
    "            dNdz = shape_gradient(fe_values, q_point, i)[3]\n",
    "\n",
    "            B[1, i * 3-2] = dNdx\n",
    "            B[2, i * 3-1] = dNdy\n",
    "            B[3, i * 3-0] = dNdz\n",
    "            B[4, 3 * i-1] = dNdz\n",
    "            B[4, 3 * i-0] = dNdy\n",
    "            B[5, 3 * i-2] = dNdz\n",
    "            B[5, 3 * i-0] = dNdx\n",
    "            B[6, 3 * i-2] = dNdy\n",
    "            B[6, 3 * i-1] = dNdx\n",
    "        end\n",
    "        \n",
    "        A_mul_B!(DB, Ee, B)\n",
    "        At_mul_B!(BDB, B, DB)\n",
    "        scale!(BDB, getdetJdV(fe_values, q_point))\n",
    "        for p in 1:size(Ke,1)\n",
    "            for q in 1:size(Ke,2)\n",
    "                Ke[p, q] += BDB[p, q]\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    return Ke\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stiffness using scalar values\n",
    "function ke_element!{T,dim}(Ke, X::Vector{Vec{dim, T}}, fe_values::CellScalarValues{dim}, C)\n",
    "    n_basefuncs = getnbasefunctions(fe_values)\n",
    "    @assert length(X) == n_basefuncs\n",
    "    reinit!(fe_values, X)\n",
    "    @inbounds for q_point in 1:getnquadpoints(fe_values)\n",
    "        for a in 1:n_basefuncs\n",
    "            for b in 1:n_basefuncs\n",
    "                ∇ϕa = shape_gradient(fe_values, q_point, a)\n",
    "                ∇ϕb = shape_gradient(fe_values, q_point, b)\n",
    "                Ke_e = dotdot(∇ϕa, C, ∇ϕb) * getdetJdV(fe_values, q_point)\n",
    "                for d1 in 1:dim, d2 in 1:dim\n",
    "                    Ke[dim*(a-1) + d1, dim*(b-1) + d2] += Ke_e[d1,d2]\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return Ke\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stiffness using vector values\n",
    "function ke_element2!{T,dim}(Ke, X::Vector{Vec{dim, T}}, fe_values::CellVectorValues{dim}, C)\n",
    "    n_basefuncs = getnbasefunctions(fe_values)\n",
    "    @assert length(X) * dim == n_basefuncs\n",
    "    reinit!(fe_values, X)\n",
    "    ɛ = [zero(SymmetricTensor{2, dim, T}) for i in 1:n_basefuncs]\n",
    "    @inbounds for q_point in 1:getnquadpoints(fe_values)\n",
    "        for i in 1:n_basefuncs\n",
    "            ɛ[i] = symmetric(shape_gradient(fe_values, q_point, i)) \n",
    "        end\n",
    "        dΩ = getdetJdV(fe_values, q_point)\n",
    "        for i in 1:n_basefuncs\n",
    "            ɛC = ɛ[i] ⊡ C\n",
    "            for j in 1:n_basefuncs\n",
    "                Ke[i, j] += (ɛC ⊡ ɛ[j]) * dΩ\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return Ke\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = 200e9\n",
    "ν = 0.3\n",
    "λ = E*ν / ((1 + ν) * (1 - 2ν))\n",
    "μ = E / (2(1 + ν))\n",
    "δ(i,j) = i == j ? 1.0 : 0.0\n",
    "g(i,j,k,l) = λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))\n",
    "\n",
    "C = SymmetricTensor{4, 3}(g)\n",
    "\n",
    "\n",
    "M = λ/ν * (1 - ν)\n",
    "\n",
    "Cmat = [ M      λ      λ    0.0    0.0   0.0;\n",
    "         λ      M      λ    0.0    0.0   0.0;\n",
    "         λ      λ      M    0.0    0.0   0.0;\n",
    "        0.0    0.0    0.0    μ     0.0   0.0;\n",
    "        0.0    0.0    0.0   0.0     μ    0.0;\n",
    "        0.0    0.0    0.0   0.0    0.0    μ]\n",
    "\n",
    "\n",
    "interpolation = Lagrange{3, RefCube, 1}()\n",
    "quad_rule = QuadratureRule{3, RefCube}(1)\n",
    "values = CellScalarValues(quad_rule, interpolation);\n",
    "vector_values = CellVectorValues(quad_rule, interpolation);\n",
    "\n",
    "# Generate some coordinates\n",
    "x = [-1.0 -1.0 -1.0;\n",
    "      1.0 -1.0 -1.0;\n",
    "      1.0  1.0 -1.0;\n",
    "     -1.0  1.0 -1.0;\n",
    "     -1.0 -1.0  1.0;\n",
    "      1.0 -1.0  1.0;\n",
    "      1.0  1.0  1.0;\n",
    "     -1.0  1.0  1.0;]\n",
    "x = x .+ 0.05 * rand()\n",
    "x_vec = reinterpret(Vec{3, Float64}, x, (8,));\n",
    "\n",
    "n_basefunctions = getnbasefunctions(vector_values)\n",
    "Ke = zeros(n_basefunctions, n_basefunctions)\n",
    "Ke2 = copy(Ke)\n",
    "Ke3 = copy(Ke)\n",
    "\n",
    "B   =  zeros(6, n_basefunctions)\n",
    "DB  =  zeros(6, n_basefunctions)\n",
    "BDB =  zeros(n_basefunctions, n_basefunctions);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fill!(Ke, 0)\n",
    "fill!(Ke2, 0)\n",
    "fill!(Ke3, 0)\n",
    "ke_element!(Ke2, x_vec, values, C)\n",
    "ke_element2!(Ke3, x_vec, vector_values, C);\n",
    "ke_element_mat!(Ke, x_vec, values, Cmat, B, DB, BDB);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[1m\u001b[32mTest Passed\n",
       "\u001b[39m\u001b[22m"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Base.Test\n",
    "@test norm(Ke - Ke2) / norm(Ke) < 1e-14\n",
    "@test norm(Ke - Ke3) / norm(Ke) < 1e-14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stiffness successful\n"
     ]
    }
   ],
   "source": [
    "println(\"Stiffness successful\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "julia-0.6",
   "display_name": "Julia 0.6.0",
   "language": "julia"
  },
  "language_info": {
   "mimetype": "application/julia",
   "file_extension": ".jl",
   "version": "0.6.0",
   "name": "julia"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
